Bosonic reaction-diffusion processes on scale-free networks 



Andrea Baronchelli, Michele Catanzaro and Romualdo Pastor-Satorras 
Departament de Fisica i Enginyeria Nuclear, Universitat Politecnica de Catalunya, Campus Nord B4, 08034 Barcelona, Spain 

(Dated: February 22, 2008) 

Reaction-diffusion processes can be adopted to model a large number of dynamics on complex 
networks, such as transport processes or epidemic outbreaks. In most cases, however, they have 
been studied from a fermionic perspective, in which each vertex can be occupied by at most one 
particle. While still useful, this approach suffers from some drawbacks, the most important probably 
being the difficulty to implement reactions involving more than two particles simultaneously. Here 
we introduce a general framework for the study of bosonic reaction-diffusion processes on complex 
networks, in which there is no restriction on the number of interacting particles that a vertex 
can host. We describe these processes theoretically by means of continuous time heterogeneous 
mean-field theory and divide them into two main classes: steady state and monotonously decaying 
processes. We analyze specific examples of both behaviors within the class of one-species process, 
comparing the results (whenever possible) with the corresponding fermionic counterparts. We find 
that the time evolution and critical properties of the particle density are independent of the fermionic 
or bosonic nature of the process, while differences exist in the functional form of the density of 
occupied vertices in a given degree class k. We implement a continuous time Monte Carlo algorithm, 
well suited for general bosonic simulations, which allow us to confirm the analytical predictions 
formulated within mean-field theory. Our results, both at the theoretical and numerical level, can 
be easily generalized to tackle more complex, multi-species, reaction-diffusion processes, and open a 
promising path for a general study and classification of this kind of dynamical systems on complex 
networks. 

PACS numbers: 89.75.-k, 87.23. Ge, 05.70.Ln 



I. INTRODUCTION 

Many natural, social, and artificial systems exhibit 
heterogeneous patterns of connections and interactions 
that can be naturally described in terms of networks or 
graphs Thus, complex network theory turns out to be 
the natural framework in which the functional and struc- 
tural properties of complex systems belonging to com- 
pletely different domains can be rationalized and inves- 
tigated [2j, y, |j, |5| . This approach has recently proved to 
be very powerful, and systematic statistical analysis have 
allowed to recognize the existence of many characteristic 
features shared by a large class of different systems, the 
most peculiar being the small-world property [f| and a 
large connectivity heterogeneity yielding a scale-free de- 
gree distribution [tJ. A graph is said to be small- world 
when the average topological distance between any pair 
of vertices is "small", scaling logarithmically or slower 
with the system size N. On the other hand, defining the 
degree k of a vertex as the number of connections linking 
it to other vertices, scale-free (SF) networks are charac- 
terized by a degree distribution P(k) that decreases as a 
power-law, 



P(k) ~ fc" 7 , 



(1) 



where 7 is a characteristic degree exponent, usually in 
the range 2 < 7 < 3. 

The distinctive structural properties of networked sys- 
tems, beyond being intrinsically interesting, have also a 
strong impact on the dynamical processes taking place 
on such systems Q, which can have practical implica- 
tions in, e.g., understanding traffic behavior in techno- 



logical systems such as the Internet [9(. In particular, 
the heterogeneous connectivity pattern of SF networks 
with diverging second moment (k 2 ) (i.e., with 7 < 3) can 
lead to very surprising dynamical properties, such as an 
extreme weakness in front of targeted attacks, aimed at 
destroying the most connected vertices fiol. Ill, as well 



as to the propagation of infective agents 12 , 13| . After 
those initial discoveries, a large series of new results has 
been put forward and we refer the reader to Refs. [1, [T3| 
for recent reviews on the subject. 

A powerful framework to describe many dynamical 
processes in a most general way is given by the the- 
ory of reaction-diffusion (RD) processes [HI]. RD pro- 
cesses are defined in terms of different kinds of parti- 
cles or "species" , which diffuse stochastically (usually by 
performing a random walk) and interact among them ac- 
cording to a given set of reaction rules. Apart from their 
natural application to describe chemical reactions, RD 
processes are useful to represent any system in which dif- 
ferent kinds of "agents" diffuse in space and dissappear, 
are created, or change their state, according to the state 
of other agents in a given neighborhood. An example of 
this kind of processes is the spread of diseases in pop- 
ulation systems. An epidemic process leading to an en- 
demic state can be described by the Susceptible- Infected- 
Susceptible (SIS) model [la ], which corresponds to an 
RD process with two species (individuals), representing 
susceptible (S) and infected (I) individuals, which dif- 
fuse (with possible different diffusion rates) and interact 
through the reactions 



S + I 
I 



2/, 
5, 



(2) 
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representing susceptible individuals becoming infected 
upon encountering an infected individual, and infected 
individuals spontaneously recovering. 

RD processes on regular topologies (Euclidean lattices) 
have been extensively studied, and an elegant formalism 
has been developed to allow for a general description in 
terms of field theories On the other hand, the ef- 

fects of more complex, heterogeneous, topologies have 
been taken into account only recently for simple pro- 
cesses QJ, QJ, HO, HI El HI, and a systematic descrip- 
tion of this interesting problem is still lacking. More- 
over, so far most of the attention has been devoted to 
the restricted case of fermionic (or microscopic, in the 
chemistry jargon) RD processes, in which a vertex of 
the network cannot be occupied by more than one par- 
ticle. In this context, numerical and analytical results 
have been put forward for the most simple RD pro- 
cesses, namely the diffusion-annihilation [IH, and the 
diffusion-coagulation [2lJ processes. Although these 
results are undoubtedly interesting and offer an initial 
insight into the behavior of RD processes in heteroge- 
neous networks, the adopted fermionic approach suffers 
from two considerable conceptual drawbacks: (i) there is 
no systematic framework for the description of this kind 
of processes, and both numerical models and theoretical 
approximations (through heterogeneous mean-field the- 
ory) must be considered on a case by cases basis; and 
(ii) it is relatively easy to deal with RD processes with 
at most order-two reactions (involving at most two par- 
ticles), but it becomes more problematic to implement 
reactions among three or more particles. Thus, for ex- 
ample, a fermionic study of the three particles reaction 
A + B + C — ► [H[ requires the introduction of an artifi- 
cial "intermediate" particle, created from the reaction of 
two particles, and that reacts itself with a third, leading 
to the actual annihilation event. In this sense, it seems 
more natural and realistic to consider instead bosonic (or 
mesoscopic) processes, in which there are no restrictions 
on the vertex occupancy, and for which, levering in what 
is already known for Euclidean lattices, it is possible to 
develop systematic analytical and numerical formalisms. 
Moreover, while some processes may naturally fit into a 
fermionic framework, other are intrinsically bosonic. For 
example, during the spreading of a disease (say HIV) 
on a social interaction network, each individual can only 
change its state (become infected) by contagion through 
one acquaintance in the network. However, the diffusion 
of a disease at the level of airport networks (for example, 
SARS) is better modeled by taking into account the num- 
ber of infected individuals in each city [24| ■ The bosonic 
version of RD processes on complex networks has been 
so far neglected, with the exception of Ref . [25| , where it 
has been applied to the particular case of the SIS process, 
Eq. (J2) (see also Ref. [26| for an extension of the model 
to weighted networks). 

In this paper, we investigate the properties of gen- 
eral bosonic RD processes in complex heterogeneous net- 
works, adopting a twofold continuous time approach 



based on heterogeneous mean-field (MF) theory and nu- 
merical simulations. We develop a general MF formalism, 
based on the standard law of mass action, that is able to 
describe any RD processes on general complex networks, 
taking a particularly simple form in one-species RD sys- 
tems. For this case, general predictions, independent of 
the particular form of the reaction rules, can be made in 
the small particle density (diffusion-limited) regime. The 
formalism is applied and fully solved in two particular 
cases, the branching-annihilating random walk and the 
diffusion-annihilation problem, examples of RD systems 
with stationary states and monotonously decaying parti- 
cle densities, respectively. In order to check the possible 
differences between the bosonic and fermionic implemen- 
tations of the same problem, we consider at the same time 
both examples from the fermionic MF theory perspective. 
We find that both formalisms provide analogous results 
for the time evolution and critical properties of the dy- 
namics. However, the two approaches are not completely 
equivalent: the functional form of the particle density 
restricted to vertices of given degree k varies widely be- 
tween the two approaches. Finally we check our results, 
and in particular the equivalence between fermionic and 
bosonic formalisms, by means of extensive c omp uter sim- 
ulations. Contrarily to previous approaches [251 ], in which 
a parallel updating scheme was defined for the particular 
model under scrutiny, we adopt a sequential continuous 
time algorithm that can be easily generalized for any RD 
process. 

The paper is organized as follows. In Sec. |TT] we de- 
fine general bosonic RD processes in complex networks. 
Sec. IIIII is devoted to the introduction of a generic ana- 
lytical framework based on bosonic heterogeneous MF 
formalism, from which general predictions can be ob- 
tained for any kind of RD process. In Sec. IIVI we con- 
sider and solve some particular examples of RD processes, 
exhibiting steady states and a monotonously decaying 
density, namely the branching annihilating random walk 
and diffusion- annihilation processes, respectively. The 
predictions of heterogeneous MF theory are validated in 
Sec. |V] by means of numerical simulations. Finally, in 
Sec. IVII we summarize and discuss our results. 



II. BOSONIC RD PROCESSES IN COMPLEX 
NETWORKS 

We consider RD processes on complex networks, which 
arc fully defined by the adjacency matrix a.y, which takes 
the values ay = 1 if vertices i and j are connected by 
an edge, and zero otherwise. From a statistical point of 
view, the network can also be described by its degree 
distribution P{k) and its degree correlations, given by 
the conditional probability P(k'\k) that a vertex of de- 
gree k is connected to a vertex of degree k' Both 
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descriptions are related through the formulas 



where 5(x, y) is the Kronecker 8 symbol, and [l| 



(4) 



TV being the size of the network. 

RD processes are defined as dynamical systems involv- 
ing particles of S different species A a , a = 1, . . . , S, 
that diffuse stochastically on the vertices of the network 
and interact among them upon contact on the same ver- 
tex, following a predefined set of R reaction rules. In a 
bosonic scheme, there is no limitation in the number of 
particles that a vertex can hold, therefore the occupa- 
tion numbers nf(t), denoting the number of particles of 
species A a in vertex i at time t, can take any value be- 
tween and oo. We will assume that diffusion in the net- 
work is homogeneous and takes place by means of random 
jumps between nearest neighbors vertices. Therefore, an 
A a particle with a diffusion coefficient D a at vertex i will 
jump with a probability per unit time D a /ki to a vertex 
j adjacent to i, where ki is the degree of the first vertex. 

The reaction rules that particles experience upon con- 
tact, on the other hand, can be defined in the most gen- 
eral way by the corresponding stoichiometric equations 

M 



s 

E 

a=l 



q r a A a ^J2(q r a +p r a )A a , r = 1 , 



, R; 



(5) 



a=l 



where q r a > (we do not consider reactions involving the 
spontaneous creation of particles) and p r a > —q 7 a . The 
coefficients q r a and p r a define the r-th reaction process, 
while A r is the probability per unit time that the reaction 
takes place. Given that the reactions take place inside 
the vertices, the only variation between a RD process 
in a complex network and a regular lattices lies in the 
diffusion step. As we will see, however, this variation 
alone can induce important differences between processes 
in these two reaction substrates. 



III. HETEROGENEOUS CONTINUOUS-TIME 
BOSONIC MEAN-FIELD FORMALISM 

A first analytical description of dynamical processes 
of complex networks can be obtained by means of het- 
erogeneous MF theory [8|. MF theory applied to net- 
works is based in the assumption that all vertices with 
the same degree share essentially the same dynamic prop- 
erties, and can therefore be consistently grouped into the 
same degree class. In the case of RD processes, and in 
order to allow for the possibility of network heterogene- 
ity and large degree fluctuations, it becomes necessary 



to work with the density spectra p a ,k(t) PJL [29|, repre- 
senting the partial density of A a particles in vertices of 
degree k, and that is defined as 



Pa,k(t) = 



n a ,k(t) 



(6) 



where n a ^{t) is the average occupation number of par- 
ticles A a in the class of vertices of degree k and Nk — 
NP(k) is the number of vertices of degree k in a network 
of size N. From the density spectra, the total density of 
A a particles is given by 



Pa(t) =E P ( fc K.fcW- 
k 



(7) 



Heterogeneous MF theory is given in terms of rate 
equations for the variation of the partial densities p a ,k(t), 
which in this case are composed by two terms: one deal- 
ing with the (linear) diffusion and another with the re- 
actions, so we can write 



dp a ,k(t) 
dt 



V a +1l a . 



(8) 



The diffusion term is easy to obtain by considering the 
diffusion dynamics at the vertex level. The total change 
of A a particles at vertex i is due to the outflow of particles 
jumping out at rate D ai plus the inflow corresponding to 
jumps of particles from nearest neighbors. Therefore, the 
diffusive component at the single vertex level satisfies the 
rate equation 



dt 



-D a n ati (t)+D a J2^n a>j (t), (9) 



Considering the density spectrum as the average 



Pa,k(t) = 



(10) 



and assuming that n a ^(t) ~ n a ^{t) 1 Vi G k, we obtain 

P(k'\k) 



V a = -D a p a ,k(t) + D a k 



k' 



- Pa ,k>(t), (If) 



where we have used Eq. (j4|). 

The reaction term can be directly derived from the law 
of mass action, according to which the rate of any (chemi- 
cal) reaction is proportional to the product of the concen- 
trations (or densities) of the reactants [3(| • Considering 
the set of all allowed processes Eq. ([SJ, we obtain: 



\K\{[p (i .k{t)] q ^ 



(12) 



Collecting all terms, the rate equations for the density 
spectra can be written in the most general case as 



dp a ,k{t) 
dt 



-D a p a dt) + £> a fcV \ - p a ,k'(t) 



k' 



p r a Kl[[p0,k(t)} q 



(13) 
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while the total densities satisfy the equations 



dp a (t) 

at 
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(14) 



where we have used the degree detailed balance condition 



kP(k)P{k'\k) = k'P(k')P(k\k'). 



(15) 



It is noteworthy that Eq. (|14p is explicitly independent 
of the particular form of the network's degree correla- 
tions, which only appear implicitly through the form of 
the density spectra p a ,k- 

In the following, we will focus in the analysis of one- 
species RD processes, in which a single class of particles 
diffuse and react in the system, i.e. S = 1. In this case, 
reactions of the same order can be grouped, and Eqs. [[Tc 
and (|14p take the simpler forms, omitting the a index, 



dt 

dpit) 
dt 



p{k'\k) 



-Pk 



,(t) + £r> fc (i)]«, (16) 



9>0 



= p(*)+£ r «E p (*)i/%W] 9 ' 

q>Q k 



where 



r, = -%,i) + £ P »-A r % r ,9), 



(17) 



(18) 



and we have absorbed the diffusion rate D into a redefi- 
nition of the time scale and the reaction rates A r . 

RD processes with non diverging solutions for Eqs. (|16[) 
and (fl~7|) can be generally grouped in two classes: those 
yielding a particle density monotonously decaying in time 
and those exhibiting one or more steady states, with 
possibly associated phase transitions between different 
steady states. We will examine more closely these two 
cases in the following subsections. While a full theoreti- 
cal analysis requires detailed information about the par- 
ticular form of the reactions involved and the network's 
degree correlations, it is possible, however, to make very 
general statements, and to obtain the asymptotic form of 
the solutions when the particle density p is very small. 



A. Steady-state Bosonic RD processes 

RD processes with steady states possess nonzero solu- 
tions for the long time limit of Eq. (|16p . In particular, 
imposing d t pk = 0, the steady states correspond to the 
solutions of the algebraic equation 



k P(k'\k) , . 
P k = -frl^—p—P*® 

k' 



q>l 



(19) 



Pk = is a solution of Eq. ([19]) . This equation is ex- 
tremely difficult to solve for a general correlation pattern 
P(k'\k), in order to find nonzero solutions. The condi- 
tion for this nonzero solution to exist, however, can be 
obtained for any correlation pattern by performing a lin- 
ear stability analysis [3l[ in Eq. (|16[). Neglecting higher 
order terms, Eq. (fT6|) becomes 



dpk(t) 
dt 



£ifcfc'Pfc'(i), 



where we have defined the Jacobian matrix 
Lkk> = Tid(k , k) H . 



(20) 



(21) 



It is easy to see that this matrix has a unique eigenvector 
t'k = k and a unique eigenvalue A = Tx + 1. Therefore, 
defining Ti =Ti + l = ^ r p r \ r 8{q r ,1), a nonzero steady 
state is only possible when T\ > 0, which translates in 
the presence of reaction processes with particle creation 
starting from a single particle. A p hase transition from 
a zero density absorbing state [32| can thus take place 
when Ti changes sign. It is worth noting that the tran- 
sition threshold takes the same form as in homogeneous 
MF theory, and it is thus independent of the network 
topology, contrary to what is found in the bosonic SIS 
model [251 ] . and similar to the case of the fermionic con- 
tact process (CP) [33j . This is due to the fact that SIS 
model is represented in terms of a two-species RD pro- 
cess, see Eq. @, in which, moreover, a conservation rule 
(total number of particles) is imposed. This conserva- 
tion rule, coupled to the diffusive nature of both species, 
is at the core of the zero threshold observed in the SIS on 
SF networks in the thermodynamic limit [25| . The con- 
tact process, on the other hand, belongs (in Euclidean 
lattices) to the same universality class as the one-species 
Schlogl RD process [34[ , hence the topology- independent 
threshold in the fermionic CP in networks can be under- 
stood in view of the general result just derived in the 
bosonic framework. 

To make further progress we restrict our attention to 
the case of uncorrelated networks, in which 13511 



P(k'\k) 



k'P(k') 



(k) ■ 

In this case, Eq. ([19[1 can be rewritten as 



Pk = 



<*>ri £[Ti 



(22) 



(23) 



Solving Eq. ([23p . we find an expression Pk(p), depending 
implicitly on the particle density. Inserting this solution 
into Eq. l[7[). we obtain a self-consistent equation for p, 



P 



= J2P(k)pk(p), 



(24) 



where we assume Ti ^ 0. Since we do not consider the 
spontaneous creation of particles from void (Tq = 0), 



to be solved in order to obtain p as a function of the RD 
parameters. 
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An approximate solution of Eq. (|23|) can be obtained 
in the limit of a very small particle density, that is, very 
close to the threshold. In this case, we can neglect the 
higher order terms in Eq. (f2"3"|) and obtain 



Pk 



(*>ri 



(25) 



which makes only sense for Ti < (i.e. < f i < 1, close 
to the phase transition). Inserting this expression into 
the self-consistent equation (|24[) yields no information. 
We must use, instead, the self-consistent relation coming 
from the steady-state condition of Eq. ([17]) . namely 



(26) 



9>1 



Inserting (|25| into Eq. (j26|) . and keeping only the term 
corresponding to the reactions of lowest order q m > 1, 
we obtain 



\(k"-)\r qm 



l/(9m-l) 



(27) 



where we have assumed T qm < 0. This solution indicates 
that, in a finite size network and for sufficiently small 
densities, all bosonic RD systems with an absorbing state 
show a critical point = 0, with an associated density 
critical exponent [3 — l/(q m — 1), coinciding again with 
the homogeneous MF solution. For SF networks with 
degree exponent 7 < q m + 1, the particle density is addi- 
tionally suppressed by a diverging factor (A:9 m )~ 1 /(9™~ 1 ) j 
signaling the presence of very strong size effects. For 
7 > 1m + 1, the particle density is size independent, and 
we recover the standard MF solution for homogeneous 
systems. 



B. Monotonously decaying Bosonic RD processes 



the limit case Ti = 0, that is, in the absence of one 
particle reactions. Then, in the limit pk — > 0, linear 
terms dominate in Eq. (|16p and we can write 



dp k (t) 
dt 



that is, the density behaves as in a pure diffusion prob- 
lem. The situation is thus the following: the time scale 
for the diffusion of the particles is much smaller than the 
time scale for two consecutive reaction events, therefore 
at any time the partial density is well approximated by 
a pure diffusion of particles [H, H3, [H[ , 



Pk(t) 



kp(t) 
(k) 



(29) 



proportional to the degree k and the total concentra- 
tion of particles, and independent of degree correlations. 
Inserting this quasi-stationary approximation back into 
Eq. ((Ill), we obtain 



dp® 
dt 



E 

<3>1 



■P q (t). 



(30) 



For small p, this equation is dominated by the reactions 
of smallest order q m . Therefore, assuming T qm < 0, we 
obtain the same decay in time as in the homogeneous MF 
theory, 



Pit) 



( gm -i)|r g j(fc^) 



-I/fan 



t-l/(<ln 



(31) 

again depressed by a size factor (k Qm ) 1 /( |?m x ) for 7 < 
q m + 1, and completely independent of the correlation 
pattern. 

If we were interested in the time behavior at interme- 
diate densities, finally, the full Eq. (flTf with the quasi- 
stationary approximation must be solved. This in general 
can only be done for uncorrelated networks. 



As we have seen in Sec. MI A\ a necessary condition 
for a RD system to have a decaying density is to have 
Ti < 0. In this case, since no steady states are present, 
the full Eq. (fT6| must be solved. One can proceed by 
using a quasi-stationary approximation [l9j, assuming 
dtPk(t) <C Pk(t), which will be correct at low densities 
if Pk(t) decays as a power law. Thus, neglecting the left- 
hand-side of Eq. (fl"6]) , we obtain again Eq. (|23|) . Solving 
it and inserting the corresponding expression of pk back 
into Eq. (jTTJ) , we have an approximate equation for p(t) 
that can give information about the long time behavior 
of the RD process. 

This procedure can be simplified when considering the 
limit of very large time and very small particle density, 
where the concentration of particles is so low that the RD 
process is driven essentially by diffusion. In this diffusion- 
limited regime, it is possible to estimate the behavior of 
the particle density, which turns out to be independent of 
the correlation pattern of the network. Let us consider 



IV. APPLICATIONS 

In this Section we will apply the bosonic MF formal- 
ism developed above to the study of two examples of one- 
species RD processes, the branching-annihilating random 
walk and the diffusion-annihilation processes, represen- 
tative of the classes of steady-state and monotonously 
decaying processes, respectively. For the sake of compar- 
ison, we will review also the predictions of corresponding 
fermionic MF theory, developed for an interacting parti- 
cle system defined to simulate the process under scrutiny. 



A. Steady-state processes: Branching-annihilating 
random walks 

On of the simplest RD processes leading to a nontriv- 
ial steady state is the generalized branching-annihilating 
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random walk (BARW), denned by the reactions [3£ 



qA 
A (p+l)i4 



(32) 



that is, particles annihilate in (/-tuples with a rate A, and 
produce a number p of offspring with rate p. Homoge- 
neous MF theory predicts a continuous phase transition 
at p c = 0, with a particle density in the active phase 



1/(9-1) 



(33) 



For the particular case q = 2, the transition belongs to 
different universality classes, according to the parity of 
the number of offsprings p [39(. If p is an odd number, 
it be long s to the universality class of directed percola- 
tion [32|, the same as the CP. On the other hand, an 
even p, for which the parity of the number of particles 
in conserved, leads to a new, and different, universality 
class. 



1. Fermionic MF theory 

When analyzing the process given by Eq. (|3"2"|) , the lim- 
itations of a fermionic approach become evident. Indeed, 
reactions involving more than two particles are difficult 
to describe in a fermionic framework, even from a concep- 
tual point of view. In fermionic models [H, [U HH , usu- 
ally diffusion and reactions are intimately linked, since 
particles jump between nearest neighbors and interact 
upon landing on an occupied vertex. Thus, when more 
than two particles are involved in a single reaction, com- 
plex schemes have to be devised to represent the process, 
schemes which, on the other hand, cannot be easily han- 
dled with standard sequential algorithms. Possible solu- 
tions could be the use of auxiliary "intermediate" parti- 
cles [23], or the design new algorithms that include the 
circumstance of different particles diffusing at the same 
time, but it is easy to figure out situations that would be 
potentially critical for such schemes (e.g. what happens 
in a pure diffusive process when one particle tries to move 
to an occupied vertex, while its starting point has been 
occupied by other particle?). On the other hand, it would 
be possible to construct such reaction schemes by involv- 
ing a particle and two or more of its nearest neighbors, in 
a reaction step independent of diffusion. Such formalism, 
although possible in principle, would be nevertheless not 
general, since the number of reacting particles would be 
limited by the connectivity of the considered vertex, and 
it would also be more cumbersome to analyze from a MF 
perspective. 

To allow for a consistent fermionic description, we will 
restrict our attention to the particular case q = 2, in 
which only binary annihilation events are allowed, and 
that can be defined as a fermionic interacting particle 
system given by the rules: 



• With probability /, a particle jumps to a randomly 
chosen nearest neighbor. 

— If the neighbor is empty, the particle fills it, 
leaving the first vertex empty. 

— If the neighbor is occupied, the two particles 
annihilate, leaving both vertices empty. 

• With probability 1 — /, the particle generates p 
offsprings. To do so: 

— p different neighbors are randomly chosen 

— A new offspring is created on every selected 
vertex, provided this is empty (if it is already 
occupied, nothing happens). 

In order to avoid problems with the offspring generation 
step, the minimum degree of the network is taken to be 
m > p. We note that this algorithm is not parity conserv- 
ing, but we do not expect this to be relevant in networks 
at MF level. 

With this implementation of the fermionic BARW in 
complex networks, we can see that the corresponding MF 
theory for the density spectrum takes the form 



dpk 

at 



= -fPk- fkp k j2^p(k'\k)p k ' 

k' 

+ fk(l- Pk )J2^P(k'\k)pk' 



+ (l-f)k(l-p k )J2^P(k'\k)p k >, (34) 

where p/k 1 is the probability that one offspring of a par- 
ticle in a vertex of degree k' arrives at a given nearest 
neighbor. For the particular case of uncorrelated net- 
works, this equation simplifies to 

-dr = - pk -jk) Pk { Pk){ ] W (35) 

where we have rescaled the time and defined v = (1 — 
f)p/f- The steady-state condition d t pk = yields the 
expression 



Pk 



k(l + u)p/(k) 
l + k(2 + u)p/{k)' 



(36) 



Application of the self-consistent condition p = 
Efe P(k)pk yields 



T P(k)k(l + V )p/(k) = 



(37) 



The condition for the existence of a nonzero solution, 
\l/'(0) < 1, yields the threshold for the existence of a 
steady state 



• Each vertex can be occupied by at most one particle 



v > v c 







f<fc 



1. 



(38) 



7 



In order to obtain the asymptotic behavior of p as a func- 
tion of v in infinite SF networks, we proceed to integrate 
Eq. (|37p in the continuous degree approximation, replac- 
ing sums by integrals and using the normalized degree 



distribution P(k) 



_1 (7 — l)k 7 , where m is the 



minimum degree in the network, to obtain 



1 



P = 



-F[l,7-1, 7 , 



(k) 



m(2 + v)p 



(39) 



where F[a, b, c, z] is the Gauss hypergeometric func- 
tion |40j |. Expanding the hypergeometric function in the 
limit of small p, close to the absorbing phase, we recover 
at lowest order for 7 > 3 the homogeneous MF result 
p ~ v. For 2 < 7 < 3, on the other hand, we obtain 



P ' 



,1/(7-2) 



(40) 



corresponding to an absorbing state transition, given by 
the control parameter with zero threshold and a criti- 
cal exponent (3 = 1/(7 — 2). 

In any finite network this behavior is modified by 
finite size effects. To analyze it, we define 6 = 
J2k kP(k)pk/ (k). The equation for the total density be- 
comes then 



dp 
dt 



p[v-(2 + v)Q\. 



(41) 



By imposing stationarity {dtp — 0) and non-zero solution 
(p 7^ 0) one obtains 



e = 



(2 + 1/)' 



(42) 



The expression of pk, Eq. (|36|) . can be simplified in the 
small density regime (p <§C (fc)/[fc(2 + v)\, Vfc) as 



Pk 



k{l + v)p 
(k) ■ 



(43) 



By substituting this expression in the definition of Q and 
inserting it into Eq. (|4"2"|) one obtains 



(fc 2 ) (l + i/)(2 + i/)' 



(44) 



SF networks of finite size have a cutoff or maximum de- 
gree k c (N) which is a function of N[35(. Therefore, for 
uncorrelated SF networks with degree cutoff scaling with 
the network size as k c (N) ~ A 1 / 2 , finite size effects in 
the fermionic BARW lead to a size dependent density 
scaling as 



P 



.- -(3— y) 
N 5 V. 



(45) 



2. Bosonic MF theory 



A bosonic formalism imposes no practical restriction 
to the maximum order that the reaction steps may 



have. Thus, the general BARW defined by the reactions 
Eq. (f3"2"| yields, within the bosonic MF formalism, to a 
rate equation Eq. (|16[) with Ti = pp and T q = —qX, and 
T q i = 0, for q' =/= {I,?}, corresponding to an absorbing 
state phase transition at a critical particle creation rate 
p c = 0. The full analysis of this equation for any q can be 
cumbersome, but we can immediately predict the behav- 
ior at large times in finite networks, which will be given 
by Eq. (|2?J). namely 



pp)Y- 



q \ 1/(9-1) 



,1/(9-1) 



(46) 



for uncorrelated networks. 

To proceed further, we consider the simplest case q 
2, in which the density spectrum fulfills the equation 

kp 



\UpI 

yielding the solution 

Pk = 



(k) 



0. 




(47) 



(48) 




where, in order to ensure the existence of the absorbing 
state, we must impose the condition Ti < 0. In the large 
kp regime, we observe here a distinctively square root 
dependence, 



(49) 



different from the limiting constant behavior observed 
in the corresponding fermionic formulation, Eq . (|36p . as 
well as in other fermionic models [12, IH, l33| . In the 
low density regime, on the other hand, we can Taylor 
expand Eq. |48p and recover, for the particular case of the 
BARW, the general relation Eq. (|25| . Thus, for particle 
densities smaller than the crossover density p x , with 



4|r 2 | /0x fc c 
(fc)|r x p 



(50) 



we recover, for uncorrelated SF networks, the asymptotic 
finite size solution for q m — 2, given by Eq. (|46|) . 

For networks in the infinite size limit, introducing the 
density spectrum of Eq. (|48|) into the self-consistent equa- 
tion Eq. (|24p. we obtain 




(51) 



In the continuous degree approximation, we have 

2(7- 



P 



2|r 2 



-1 



27- 



'4|T2|mp 

Wiril 2 



[ 2' 7 2' 7 2' A\T 2 \mp l 



(52) 
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Expanding F[a, b, c, z] in the limit of small p, we find 



• Each vertex can be occupied by at most one particle 



iril 



|ri| 4|r 2 |V^ 
l(fc)(ri) 2 



r(2- 7 )r( 7 -3/2)x 

0{ P 2 ). (53) 



At lowest order, and for 7 > 3, we recover the homoge- 
neous MF solution p ~ ti ~ p/i. On the other hand, for 
2 < 7 < 3, the nonzero solution of this equation is 



(pa*) 



1/(7-2) 



A 



(54) 



corresponding to an absorbing state transition, given by 
the control parameter p, with zero threshold and a crit- 
ical exponent [3 = 1/(7 — 2), in full agreement with 
the results for the corresponding fermionic version of 
the model. We can use this last result to estimate the 
crossover density to the finite size solution Eq. (|46|) . In- 
serting Eq. ((54} into Eq. ([50]) . and considering T2 as a 
constant, we obtain that the finite size solution should 
be observed for a control parameter 



p < p> 



V 



(55) 



Therefore, for uncorrelated SF networks, finite size ef- 
fects in the bosonic BARW should appear for a particle 
creation rate smaller that p x ~ jV^ 7-2 )/ 2 . 



B. Decay processes: Diffusion-annihilation process 

The simplest case in the class of monotonously decay- 
ing RD processes corresponds to the general diffusion- 
annihilation process 



qA 



(56) 



which is the particular case of the BARW analyzed in 
Sec. IIV Al with p, = (at the critical point). The ho- 
mogeneous MF solution predicts a decay of the particle 
density 



^(i)~t- 1/ <«- 1 '. 



(57) 



In Euclidean lattices of dimension d, dynamical renor- 
malization group arguments [41| show that the behavior 
in Eq. (|5T[) is correct for d above the critical dimension 
d c = 2/(q — 1). Below it, we have instead p(t) ~ t~ d / 2 , 
with logarithmic corrections appearing at d — d c . 



Fermionic MF Theory 



As discussed in Sec. IIV A 11 in order to allow for a 
consistent fermionic description, we will restrict our at- 
tention to the binary diffusion-annihilation process with 
<7 = 2, which can be implemented as a fermionic inter- 
acting system obeying the following rules [ill : 



• Each particle jumps with probability / to a ran- 
domly chosen nearest neighbor. 

• If the neighbor is empty, the particle fills it, leaving 
the first vertex empty. 

• If the neighbor is occupied, the two particles anni- 
hilate, leaving both vertices empty. 

This model was analyzed in detail in Ref. [Hj]. There 
it was observed that the rate equation for the density 
spectrum reads, in uncorrelated complex networks, 



dpk 
dt 



= -Pk 



(58) 



where the probability / has been absorbed into a rescal- 
ing of time. With a quasi-stationary approximation, the 
density spectrum at large times takes the form 



Pk{t) 



kp(t)/(k) 



(59) 



1 + 2kp(t)/(k) ' 
which yields as a final equation for the density of particles 



dp 
di 



» 



k 



P{k) 



1 + 2kp(t)/(k) ' 



(60) 



In finite networks, for times larger that t > i x , with 
k c p(t x ) ~ 1, the denominator in Eq. (|60|) can be sim- 
plified to 1, to obtain the limit behavior in finite size 
networks 



P(t) 



2(k 2 



r 



(61) 



In a network of infinite size, the full Eq. (|6"0")) must be 
integrated. Within the continuous degree approximation, 
this equation takes the form 

^ = -p(t)F[l, 7 - 2, 7 - l,-(k)/2mp(t)}. (62) 

Expanding the Gauss hypergeometric function for small 
p, we obtain, for 7 > 3, the asymptotic long time behav- 
ior pit) ~ t^ 1 while for 2 < 7 < 3, one has 



p{t)~t- 1 ^-*). 



2. Bosonic MF Theory 



(63) 



The general diffusion-annihilation process defined by 
reaction Eq. (|56p leads to the general rate equation 
Eq. p6[) . with parameters T\ = 0, T q = —q\, and 
Tg> = 0, for q' 7^ {l,?}- In finite networks and for large 
times, the behavior of the particle density will be given 
by Eq. (ST]), i.e. 



Pit) 



(q-l)qX(k") 



-1/(9-1) 



-1/(9-1) 



i-V(9-i) 



(64) 
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the last expression holding for uncorrelated SF networks. 

Let us focus again in the simplest case q = 2. The rate 
equation for the total density in uncorrelated networks 
takes the form 



dp(t) 

at 



\T 2 \J2P(k)p 2 k (t). 



(65) 



Applying the quasi-stationary approximation for the den- 
sity spectrum, we are led to the second order equation 



\F2\pl+Pk- 7TT P = Q 



(*>' 



whose only positive solution is 



Pk 



2 r 



1 '-i + 4 / 1+ 4£!V 



(66) 



(67) 



For a finite network with degree cut-off fc c , when the 
density is smaller that 



(68) 



we can Taylor expand Eq. (|67|) to obtain the expres- 
sion pk ~ kp/(k) and the asymptotic behavior given by 
Eq. (|64|) . On the other hand, for large k and p, we obtain 



Pk 



kp 



(69) 



and we find again the peculiar square root behavior of 
the density spectrum on k, distinctive from the fermionic 
prediction. 

The general solution in the infinite network limit 
can be obtained in this case by substituting the quasi- 
stationary approximation (|67[) into Eq. (|65[) . to obtain 



i-(i + 2|r 2 | P ) + |r 2 |£p(fe)Ji 



4|r 2 | fc 

(k) 



-/>■ 



dp 
dt = 

(70) 

In the continuous degree approximation, and for SF net- 
works, we obtain in the infinite network size limit 



dp 



1 



7-1 U\T 2 \mp 



dt ~ 9 2|r 2 | |r 2 |(2 7 -3)y (A) 
1 3 1 ( k ) 1 

x *T-2>7- 2-7- 2'"4|r 2 |m/9 J ' 



(71) 



Considering the limit of large times and small densities, 
we can expand the hypergeometric function [401 ] . to ob- 
tain 



dp r(2- 7 )r( 7 -3/2) {4m\T 2 \p 



7-1 



-CIV) (72) 



dt ~ 4|r 2 |V^ V (*) 

for 2 < 7 < 3, whose solution is 

p(t) ~ \Y 2 \-<~H- 1 /^ ~ A 7_2 t~ 1 ^ 7_2) (73) 



that is, a power law decay with an exponent l/( 7 — 2), 
again in agreement with the fermionic implementation 
of the process. From this expression we can estimate the 
time at which the crossover density in Eq. (|6"8|) is reached 
in SF network, namely 



t x ~ fcr 2 ~ iv (7 - 2)/2 , 



(74) 



taking the same functional form as the crossover control 
parameter for BARW in Eq. (|55|) . 



V. NUMERICAL SIMULATIONS 

As we have seen in the previous Sections, het- 
erogeneous MF theory applied to steady state and 
monotonously decaying bosonic RD processes can make 
general predictions for the asymptotic behavior at finite 
networks, as well as give specific solutions for the in- 
finite network size limit. In particular, we have seen 
that bosonic formalisms provide exactly the same results 
as their fermionic counterpart (whenever the fermionic 
mapping is possible) regarding the evolution of the par- 
ticle density, the only difference being the form of the 
density spectra as a function of the degree k. In order to 
check these conclusions, we have performed extensive nu- 
merical simulations of bosonic and fermionic versions of 
the processes considered. To generate the network sub- 
strate for the RD processes, we have adopted the uncorre- 
lated configuration model (UCM) [42] that has the double 
benefit of producing SF networks without degree correla- 
tions [IE EH and with a tunable degree exponent. When 
correlations were desired, the configuration model (CM) 
HE EE EE El] was used, with the additional constraint 
of lack of multiple connections and self- loops [IE EE HE] • 

Numerical simulations of fermionic RD processes must 
be tailored on a case by case basis, depending on the 
specific interacting particle system chosen to represent 
it [IE EE EH- As a general rule, simulations are per- 
formed following a sequential Monte-Carlo scheme [32j |. 
At the beginning, Npo particles are randomly distributed 
on the network, respecting the fermionic constrain that 
at most one particle can be present on a single vertex, 
Le. po < 1. Then, at time t, a particle is randomly 
selected, and it undergoes the corresponding stochastic 
dynamics. The system is then updated according to the 
actions performed by the selected particle, and finally 
time is increased as t — ► t + l/n(t), where n(t) is the 
number of particles at the beginning of the simulation 
step. For bosonic processes, we have used a continuous 
time formalism, details of which are given in the following 
subsection. 



A. Continuous time bosonic simulations 

Previous approaches to the numerical simulation of 
bosonic RD processes on complex networks [25j j relied 
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on a parallel updating rule in which reaction and diffu- 
sion steps alternate: after all vertices have been updated 
for reaction, particles diffuse. This approach, while fea- 
sible, must again be tailored in a case by case basis, and 
strongly depends on the specific reactions of the process 
under consideration. Moreover, it introduces a subtle 
but relevant problem as far as the density spectrum is 
concerned. Indeed, while preserving the average density, 
pure diffusion immediately sets up the characteristic lin- 
ear behavior Pk(t) ~ k. Thus, the density spectrum may 
assume (very) different aspects if we look at it after the 
reaction step or after the diffusion one. In order to over- 
come these difficulties we have opted instead for a se- 
quential algorithm, which not only is absolutely general, 
but is in addition closer to the spirit of the continuous 
time rate equations we have developed to describe het- 
erogeneous MF theory. The algorithm implemented is 
based in the one proposed in Refs. [FJ, [HJ for the case 
of regular lattices. For one-species RD processes, the al- 
gorithm is described as follows: In networks of size N, 
initial conditions for simulations are usually a number 
PqN of particles randomly distributed on the network 
vertices, with no limitation on the occupation number of 
single vertices. To perform the dynamics, we consider 
the microscopic configuration {C} of the bosonic system, 
which is specified by the occupation number m at each 
vertex i. A standard master equation approach [l5T ] im- 
plies that, for RD processes described by Eq. ([8]), the 
average number of events in an infinitesimal time dt is 

E(dt, {C}) =dtJ2 (V ! A r + J£IL^ ^ ef ) 

(75) 

where 



Lo(m,q r ) = 




is the number of non-ordered q r tuples of particles at ver- 
tex i. Since the algorithm considers all reacting q— tuples 
as equivalent, it is convenient focusing on reaction or- 
ders q rather than on specific reactions r. In general, 
a particular RD process defines a finite set Q of al- 
lowed reaction orders q, that can be formally indicated 
as Q = ({q} \ 3r : q r = q). At each time step one has to: 
(i) select a vertex i (ii) select the order q of the candi- 
date reaction (iii) determine which reaction r occurs. In 
details: 

(i) A vertex i is selected with probability Wi/M, where 

Wi = J2 q eQ w ("i> ?) and M = Ei Wf, 

(ii) A particular q — q* (with q* £ Q) is selected with 

probability uj(rii, q*)/Wi\ 

(iii) A particular reaction r of order q r = q* occurs with 
probability q r \X r At, where At is a configuration- 
independent time constant. In case q* = 1, in ad- 
dition to reaction processes, the particle has the 




FIG. 1: Average density of the bosonic BARW with q = p — 
2 at the steady state on uncorrelated UCM networks with 
7 = 2.5. The annihilation rate is kept fixed at A = 0.1. 
Top: Density at the stationary state as a function of fj,, for 
different network sizes N. At any value of fj,, larger network 
sizes corresponds to smaller densities. Bottom: Check of the 
collapse predicted by Eq. (|46)l . The dashed line has slope 1. 

diffusion option, which is chosen with probability 
At (since we set the diffusion coefficient D = 1). 

Time is updated as t — » t + At/M . It is clear [5l[ that 
to have valid transition probabilities At must be chosen 
so that the condition 

ls( q) l) + q\ £ A r jAi<l (77) 

\ r:q r =q J 

holds for all values of q. With this prescription an average 
of E(At, {C}) events occur in a time interval At. 

B. Branching-annihilating random walks 

In our numerical study of the BARW, we first focus 
in the behavior of the average particle density in the 
steady state as a function of the branching rate. As al- 
ready observed in other dynamical systems in SF net- 
works [III [53| , we find it difficult to observe the infinite 
size limit behavior [Eq. ([^0]) or (|54[) ] in either bosonic 
of fermionic simulations, for the network sizes available 
within our computer resources. Therefore, we report the 
results for the finite size behavior, expected in finite net- 
works, Eqs. and (|46[) . In Fig. [1] we plot the average 
density in the active phase of the bosonic BARW with 
q = p = 2 as a function of the branching rate [i. In 
the parameter range shown in this Figure (top panel), 
we observe that the density follows a linear behavior as 
a function of the branching parameter fi. This linear de- 
pendence on fi corresponds to the asymptotic finite size 
solution predicted by Eq. PS|) , which is expected to hold 
in networks of finite size and for very small steady state 
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FIG. 2: Average density of the fermionic BARW with q = 
p — 2 at the steady state on uncorrelated UCM networks with 
7 = 2.5. Top: Density at the stationary state as a function 
of the parameter v, for different network sizes N. Bottom: 
Check of the collapse predicted by Eq. (|45p . The dashed line 
has slope 1. 



densities. We can further check the accuracy of the pre- 
diction by noticing that, in SF uncorrelated networks, 
the prefactor in p should scale with the system size as 
p ~ pN~( 3 ^/ 2 . Therefore, we should expect that a plot 
of N^-^/ 2 p as a function of p would collapse for differ- 
ent network sizes. This is actually what we observe in 
Fig. [T] (bottom panel), where different curves are clearly 
laid one on top of the other for small values of p. As 
the density becomes larger, on the other hand, the col- 
lapse becomes less and less precise, in agreement with the 
fact that Eq. ([41))) is only valid in the very small density 
regime. Moreover, deviations from the collapse line set in 
earlier for large system sizes in agreement with Eq. (|55|) , 
according to which finite size effects show up for values 
of p smaller than p x ~ aMt-2)/2_ 

In Fig. [5] we present analogous results for the fermionic 
version of the BARW. Here (top panel) we can observe 
a first difference with respect to the bosonic BARW: For 
small values of N, it is not possible to span a range of 
small values of v, due to the fact the the system falls 
quickly into the absorbing state. Small v can only be 
explored using large N. The trend of all plots is, however, 
correct: Linear in v and decreasing when increasing the 
network size. The data again collapses with the same 
functional form, now p ~ i>N~( 3 ~' 1 ^ 2 , for large systems 
sizes. The deviations at small N and large v, however, 
seem now larger than in the bosonic case. 

Having checked that the average density takes the same 
form in both bosonic and fermionic approaches, we focus 
now in the density spectra, in which differences between 
the two formalisms are predicted at the MF level. In the 
case of the bosonic BARW with q = 2, the density spectra 
in the steady state, as given by Eq. P5|) , is characterized 
by a peculiar square root behavior. To check this form, 



FIG. 3: Density spectra in the bosonic BARW process with 
q = p — 2 at the steady state on uncorrelated UCM networks 
with 7 = 2.5. Network size N = 10 6 . Top: Density spec- 
tra as a function of the degree k for different steady state 
densities. Different stationary densities have been obtained 
fixing the annihilation parameter A = 0.05, and varying the 
branching parameter p. Center: Data collapse of the den- 
sity spectra with different average stationary densities as pre- 
dicted by Eq. (|78[1 . Bottom: Check of the Taylor expansion 
of the density spectra, as given by Eq. (|79|l . 



we observe that, if we define the function 



l-2p 



(l-2/i) 2 (fc) 
8Ap 



(78) 



we expect G fl (pk) = k for any values of the reaction pa- 
rameters. In Fig. [3K center panel) we can see that this 
collapse works well for a wide range of p values. Alterna- 
tively, we can consider the small density behavior, given 
by the general Eq. (f23|) , which translates in the function 



TM = (l-2p)(k) 



Pk 



(79) 



being T^(pk) = k. In Fig. [3]Jbottom panel) we observe 
a poor collapse of the curves, which is approximately at- 
tained only at very low densities, confirming the presence 
of strong nonlinearities at large p. 

In Fig. 2] we investigate the density spectrum of a 
fermionic BARW for different values of the total density. 
As we can observe (top panel), the spectra saturates to 
a constant value for large values of p and k, as expected 
from the theoretical expression Eq. ([55]) . On the other 
hand, this equation implies that the function 



G v {pk) 



(k)p k 



p{t)[(l + v)-(2 + v)p k ] 



(80) 



should satisfy G v (pk) — k for all / and p. Considering the 
small density limit, on the other hand, a linear behavior 
of pk with k is expected, translated again in the new 
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FIG. 4: Density spectra for the fermionic BARW process with 
q = p — 2 at the steady state on uncorrelated UCM networks 
with 7 = 2.5. Network size iV = 10 6 . Top: Density spectra as 
a function of the degree k for different steady state densities. 
Different stationary densities have been obtained varying the 
parameter v. Center: Data collapse of the density spectra 
for different steady state densities, as predicted by Eq. (|80|l . 
Bottom: Check of the Taylor expansion of the density spectra, 
as given by Eq. (|81|) . 



function 



T v {p k ) 



_ (k}pk 



P(*)(l 



(81) 



being T v (p k ) — k. While the collapse with the full shape 
of Eq. l[36|) (center panels in Fig. |4|) is almost perfect, 
it is much worse if only the Tailor expansion in consid- 
ered (bottom panel), being only approximately correct 
for very small densities. 



C. Diffusion-annihilation process 

To validate our theoretical approach for decaying RD 
systems, we have concentrated on the bosonic description 
of the processes, since the fermionic version described in 
Sec. lIVBll was already checked numerically in Ref. [lj|. 
We consider thus the general bosonic process qA — > 0, at 
rate A, for which a detailed analytical solution was given 
in Sec. IIVB 21 For the case q — 2, again a peculiar square 
root behavior for the density spectrum was predicted in 
Eq. (frJT)) . which is corroborated in Fig. [5] by means of 
three different graphs. Again, from Eq. (|67|) . defining 
the function 

G (p k ) ee ((4X Pk + I) 2 - 1) @L, (82) 

where A is the annihilation parameter, we will expect 
that Go(pk) — k for all times. In Fig. [5] (center panel) 
it is clear that the different curves, corresponding to dif- 
ferent values of the average density p(t), collapse well in 



FIG. 5: Density spectra of the bosonic RD process 2A — > 
on uncorrelated UCM networks with 7 = 2.5. Network size 
TV = 10 6 . Top: Density spectra as a function of the degree k 
at different times (densities) from measures performed with 
fixed parameter A = 0.1. The curves show a bending in the 
large k region for short times (large densities). Center: Data 
collapse of the density spectra at different times as predicted 
by Eq. (|82p. Bottom: Check of the Taylor expansion of the 
density spectra, as given by Eq. (|83[) . The poor collapse at 
large k confirms the presence of strong nonlinear terms at 
short times. 



agreement with the theoretical prediction. In the bot- 
tom panel, we check the general asymptotic expression 
for large times, Eq. (l2~5l) . In this case, for small values of 
8\kp/(k), we should expect the function 



T (pk) 



kp(t) 
(k) 



(83) 



to be To(/?fc) = k, which holds when the times are large 
enough, but shows a clear bending at large degrees and 
large densities, signature again of the fact that it is fun- 
damental to take into account the non-linearity of the 
spectrum. 

As in the case of the fermionic diffusion-annihilation 
process [l9l |. it turns out that the asymptotic expression 
for infinite networks of the total particle density, Eq. (|T3[) , 
is very difficult to observe numerically, due to the very 
small range of the extension of the power-law behavior. 
We have therefore focused again on the general prediction 
for finite networks, Eq. ([64]) . according to which the RD 
process qA — > shows a decay of the average density at 
large times of the form p{t) ~ t^Vt?- 1 ^ independently 
of the presence or absence of degree correlations. We 
present in Fig. [6] simulation results for three values of q, 
namely q — 2,3,4, on uncorrelated networks SF gener- 
ated with the UCM algorithm (main plot), and correlated 
SF networks generated with the CM prescription (inset). 
It is clear that the theoretical predictions are in perfect 
agreement with numerical data. This result is particu- 
larly relevant since, for q > 2, it concerns purely bosonic 
processes, which do not have a fermionic counterpart. 
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t (l+q-Y)/2 



FIG. 6: Density decay of the bosonic qA — > diffusion- 
annihilation processes in finite correlated and uncorrelated 
networks, for different values q — 2 (full lines), q = 3 (dashed 
lines) and q — 4 (dot-dashed lines). For all values of q, 
the graphs show a tail of the form t~ 1 '^ q ~ 1 \ as predicted 
in Eq. (|64p . both for uncorrelated (UCM) networks (main fig- 
ure) and correlated (CM) networks (inset). Data obtained 
from networks of size N = 10 5 with degree exponent 7 = 2.5. 
For all plots, the annihilation parameter was fixed at A = 0.04. 



FIG. 7: System size dependence of the density prefactor 
in the diffusion-annihilation process 3^4 — > 0. According 
to Eq. (|84p . we plot the prefactor A(N, 7) as a function of 
j\K 1+ 9-7)/2 g 00C L ij near behavior confirms the predic- 

tions of bosonic heterogeneous MF theory in the diffusion- 
limited regime. Data obtained from networks with degree 
exponent 7 = 2.5. For all plots, the annihilation parameter 
was fixed at A = 0.04. 



The time independent prefactor of Eq. (|64[) . moreover, 
states that the average density should be suppressed by 
the size term ((k q )/(k))~ 1 ^ q ~ 1 \ More precisely, we can 
rewrite Eq. §7} as p(t) ~ A(N,j)~ 1 /^~ 1 h- 1 ^i~ 1 \ with 

A(JV, 7 ) ~ Ar(i+9-f)/2 (84) 

in UCM networks, with cutoff k c (N) ~ N 1 / 2 . We have 
estimated the A(N, 7) values by linear fits of the p(i) _1 vs 
^— !/ (g— lj curves for the 3 A — ► process taking place on 
networks of different sizes (data not shown), and for two 
values of the degree exponent 7. We report the results in 
Fig. [3 where the scaling relation predicted by Eq. ([84]) 
is found to be in very good agreement with simulation 
data. 



VI. DISCUSSION AND CONCLUSIONS 

In this paper we have studied bosonic RD processes 
in SF networks introducing a general continuous-time 
framework which is well suited for both MF analytic cal- 
culations and computer simulations. At the MF level, 
we have developed the rate equations that character- 
ize any generic RD process. We have considered in 
particular one-species RD processes, for which MF the- 
ory provides a natural way to classify all possible RD 
schemes. We have analyzed in detail both steady state 
and monotonously decaying processes from a general 
perspective, focusing also on specific examples, namely 
the BARW and diffusion-annihilation processes. For 
processes characterized by reactions not involving more 



than two particles, we have compared the results with a 
fermionic version of the same problems, implemented in 
terms of discrete interacting particle systems. 

Beyond the obvious difference concerning the fact that 
the average density is bounded in fermionic processes 
while it is not in their bosonic version, both bosonic and 
fermionic MF formalisms render equivalent results for the 
density of particles in single-species RD processes, the 
main difference between both formalism lying in func- 
tional form of the density spectrum of particles. For high 
densities, the spectrum in bosonic systems goes in general 
as the power fc 1 /" 1 , where qM is the highest order of the 
reactions defining the RD system, while for fermionic sys- 
tems the behavior of the spectrum is in general algebraic. 
Thus, in the bosonic scheme, hubs become relatively less 
and less populated as more many-particle reactions are 
present, provided the average density is sufficiently high. 
In the very low density regime, on the other hand, the 
bosonic approach predicts spectra linear with k, a fact 
that allows to make general predictions for the behavior 
of any RD process in finite networks, which turns out to 
coincide with the homogeneous MF result, with a net- 
work size correction. This results confirms the relevance 
of finite size effects in dynamics on SF networks, already 
reported for fermionic systems [lj| [53[ , since the border 
between "high" and "low" densities is in general deter- 
mined by the pk c {N) product. 

Another interesting result concerns one-species bosonic 
RD processes with an absorbing state phase transition, 
where the critical point does not depend on the possible 
heterogeneity of the network, but is in general located 
at Ti > 0. This condition translates in the presence of 
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reaction processes with particle creation starting from a 
single particle and corresponds to the threshold indepen- 
dent of the network topology found in other fermionic 
systems [Hj]. In order to observe effects of the connec- 
tivity heterogeneity in the threshold, more complex RD 
schemes, such as those involving two or more species, 
must be considered 25]. On the other hand, the bosonic 
point of view allows to shed a different light on the value 
7 = 3 usually associated to a frontier between regular 
(7 > 3) and complex (7 < 3) behavior for dynamical 
systems on SF networks. We can readily see that the 
value 7 = 3 emerges simply from considering dynamical 
processes involving at most two particle interactions. For 
general interactions involving q particles, one will expect 
instead to obtain unusual results for 7 < q + 1 (see e.g. 

Eq. mno. 

The continuous time theoretical and numerical for- 
malisms presented for bosonic processes have been de- 



veloped in depth for the particular case of one-species 
processes, but they can be easily generalized to many 
species systems, opening thus the path to the study a 
large variety of processes of large relevance in the under- 
standing of the topological effects of complex networks 
on dynamic and transport phenomena. 
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